パッケージとデータの読み込み

必要なパッケージとデータを読み込みます。差分の差分法を実装事態はR内蔵関数であるlm()でできますが、データハンドリングと可視化のために{tidyverse}を、記述統計の出力のために{summarytools}パッケージを読み込みます。また、クラスター標準誤差（cluster standard error）計算と出力のためには、{estimatr}が必要です。インストール済みでない方はinstall.packages()やpacman::p_load()でインストールしておきましょう1。

必要なデータは講義用Slackで配布しているDay3_Data5.csvです。Harvard Dataverseからもダウンロード可能ですが、こちらはクリーニングが必要です。本講義はRの使い方に関する講義ではないため、ここではクリーニング済みのでデータを使用します。

# library(tidyverse) 
# library(summarytools)
# library(estimatr)
pacman::p_load(tidyverse, summarytools, estimatr)

# Load dataset
Bacon_df <- read_csv("Data/Day3_Data5.csv")

記述統計量の確認と変数のクリーニング

まずは、データのサイズを確認します。

dim(Bacon_df)
## [1] 8059   15

データの大きさは8059行15列です。それでは15個の変数名をnames()関数を使って出力します。これはデータフレームの列名を出力する関数です。

names(Bacon_df)
##  [1] "Assembly" "Year"     "Treat"    "FG"       "y_0"      "y_m1"    
##  [7] "y_m2"     "y_m3"     "y_m4"     "y_m5"     "y_p1"     "y_p2"    
## [13] "y_p3"     "y_p4"     "y_p5"

それぞれの変数の詳細はFouirnaies and Mutlu-Eren (2015)2を参照して頂きたいところですが、こちらでも簡単に説明しておきます。

変数名 説明
Assembly 地方議会名
Year 年
Treat 地方議会の多数派と政府与党の一致有無
FG t期の一人あたり定式補助金（対数）
y_0 t期の一人あたり特定補助金額（対数）
y_m* t - *期の一人あたり特定補助金額（対数）
y_p* t + *期の一人あたり特定補助金額（対数）

データの記述統計量を出力します。

# 記述統計量を確認
view(dfSummary(Bacon_df))
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 Assembly [character] 1. DURHAM 2. ADUR 3. ALLERDALE 4. AMBER VALLEY 5. ARUN 6. ASHFIELD 7. ASHFORD 8. AYLESBURY VALE 9. BABERGH 10. BARKING & DAGENHAM [ 411 others ]
38(0.5%)
21(0.3%)
21(0.3%)
21(0.3%)
21(0.3%)
21(0.3%)
21(0.3%)
21(0.3%)
21(0.3%)
21(0.3%)
7832(97.2%)
 8059 (100.0%) 0 (0.0%)
2 Year [numeric] Mean (sd) : 2001.8 (6) min < med < max: 1992 < 2002 < 2012 IQR (CV) : 10 (0) 21 distinct values 8059 (100.0%) 0 (0.0%)
3 Treat [numeric] Min : 0 Mean : 0.3 Max : 1
0:5618(72.4%)
1:2140(27.6%)
 7758 (96.3%) 301 (3.7%)
4 FG [numeric] Mean (sd) : 16.9 (1.6) min < med < max: 13.9 < 16 < 20.9 IQR (CV) : 3 (0.1) 6936 distinct values 8059 (100.0%) 0 (0.0%)
5 y_0 [numeric] Mean (sd) : 2.7 (2.1) min < med < max: -2.4 < 2.3 < 7.7 IQR (CV) : 3.2 (0.8) 7813 distinct values 7993 (99.2%) 66 (0.8%)
6 y_m1 [numeric] Mean (sd) : 2.7 (2) min < med < max: -2.4 < 2.2 < 7.7 IQR (CV) : 3.2 (0.8) 7362 distinct values 7533 (93.5%) 526 (6.5%)
7 y_m2 [numeric] Mean (sd) : 2.6 (2) min < med < max: -2.4 < 2.1 < 7.7 IQR (CV) : 3.1 (0.8) 6909 distinct values 7075 (87.8%) 984 (12.2%)
8 y_m3 [numeric] Mean (sd) : 2.5 (2) min < med < max: -2.4 < 2 < 7.4 IQR (CV) : 3.1 (0.8) 6454 distinct values 6615 (82.1%) 1444 (17.9%)
9 y_m4 [numeric] Mean (sd) : 2.4 (1.9) min < med < max: -2.4 < 1.9 < 7.3 IQR (CV) : 3.1 (0.8) 6006 distinct values 6160 (76.4%) 1899 (23.6%)
10 y_m5 [numeric] Mean (sd) : 2.4 (1.9) min < med < max: -2.4 < 1.7 < 7.3 IQR (CV) : 3 (0.8) 5584 distinct values 5735 (71.2%) 2324 (28.8%)
11 y_p1 [numeric] Mean (sd) : 2.8 (2.1) min < med < max: -2.4 < 2.3 < 7.7 IQR (CV) : 3.3 (0.7) 7397 distinct values 7561 (93.8%) 498 (6.2%)
12 y_p2 [numeric] Mean (sd) : 2.9 (2.1) min < med < max: -2.4 < 2.3 < 7.7 IQR (CV) : 3.4 (0.7) 6962 distinct values 7108 (88.2%) 951 (11.8%)
13 y_p3 [numeric] Mean (sd) : 2.9 (2.1) min < med < max: -2.4 < 2.4 < 7.7 IQR (CV) : 3.5 (0.7) 6522 distinct values 6649 (82.5%) 1410 (17.5%)
14 y_p4 [numeric] Mean (sd) : 3 (2.1) min < med < max: -1 < 2.4 < 7.7 IQR (CV) : 3.6 (0.7) 6104 distinct values 6209 (77.0%) 1850 (23.0%)
15 y_p5 [numeric] Mean (sd) : 3.1 (2.1) min < med < max: -1 < 2.4 < 7.7 IQR (CV) : 3.6 (0.7) 5707 distinct values 5783 (71.8%) 2276 (28.2%)

Generated by summarytools 0.9.9 (R version 4.0.4)
2021-08-04

Year変数は連続変数ではなく、ダミー変数として利用する予定なので、integer型からfactor型へ変換します。factor型変数を回帰モデルに投入すると自動的にダミー化して推定をしてくれるので便利です。また、後ほどトレンド変数も使いますのでトレンド変数も作成します。別にこのままでも問題ありませんが、1992年を1、1993年を2、…にしたいので、Year変数から1991を引けばいいです。

Yearをfactor型にしたら引き算ができなくなるので、トレンド変数を作成し、続いてYear変数をfactor型にします。

Bacon_df <- Bacon_df %>%
    mutate(Trend = Year - 1991,
           Year  = factor(Year))

線形トレンドを無視したDID推定

まずは、線形トレンドを考慮しない差分の差分法を行います。やり方としては、説明変数として、処置変数とユニット3のダミー変数、年度のダミー変数を投入します。記述統計で見られるように、AssemblyとYearはfactor型になっているため、このまま入れてたらダミー変数化した上で分析をやってくれます。

結果変数 (応答変数)は、\(t+1\)年の当該自治体への特定補助金 (一人あたり; 対数)から定式補助金です。ただし、同じ党派性を持つという処置から起因する補助金の増額が\(t+1\)年から見られるとは限りません。したがって、\(t+2\)から\(t+5\)年までの補助金も説明変数として使いましょう。

まずは特定補助金（y_で始まる変数）から定式補助金(FG)を引きましょう。

Bacon_df <- Bacon_df %>%
    mutate(y_0  = y_0  - FG, # y_0からy_m5は後で使います。
           y_m1 = y_m1 - FG, # 今のうちに引いておきます。
           y_m2 = y_m2 - FG,
           y_m3 = y_m3 - FG,
           y_m4 = y_m4 - FG,
           y_m5 = y_m5 - FG,
           y_p1 = y_p1 - FG,
           y_p2 = y_p2 - FG,
           y_p3 = y_p3 - FG,
           y_p4 = y_p4 - FG,
           y_p5 = y_p5 - FG)

以上のコードはacross()関数を使うと以下のように簡潔にすることができます。大変便利な関数ですので、ぜひ身につけておきたいですね。詳しい使い方は『私たちのR』の第14章を参照してください。

Bacon_df <- Bacon_df %>%
  # "y_"で始まる全ての変数からFGを引き、上書きする
  mutate(across(starts_with("y_"), ~.x - FG))

それではy_p1からy_p5までを結果変数とした差分の差分推定を行います。y_0を結果変数として使わない理由は単純です。議会と政府の党派性が一致した時点での補助金額は既に決まっているからです。補助金額に変動があるとすれば、次の年からだと考えた方が自然でしょう。

推定に投入する変数は処置変数に加え、議会と年とダミー変数となります。また、差分の差分法で使う標準誤差はクラスター標準誤差を使うのが一般的です。lm()関数にはクラスター標準誤差を計算する機能が付いておりませんので、ここでは{estimatr}のlm_robust()関数を使います4。cluster引数にはクラスタリングを行う変数名を入力します。se_typeはクラスター標準誤差の種類を指定します。選択できるオプションとしましては"CR0"、"CR2"（既定値）、"stata"のいずれかを指定できますが、ここではStataと同じアルゴリズムを使いますので、"stata"と指定します。

# Model 1: Diff-in-Diff Estimation without linear trends
Bacon_Fit1_P1 <- lm_robust(y_p1 ~ Treat + Assembly + Year,
                    data = Bacon_df, cluster = Assembly, se_type = "stata")
Bacon_Fit1_P2 <- lm_robust(y_p2 ~ Treat + Assembly + Year,
                    data = Bacon_df, cluster = Assembly, se_type = "stata")
Bacon_Fit1_P3 <- lm_robust(y_p3 ~ Treat + Assembly + Year,
                    data = Bacon_df, cluster = Assembly, se_type = "stata")
Bacon_Fit1_P4 <- lm_robust(y_p4 ~ Treat + Assembly + Year,
                    data = Bacon_df, cluster = Assembly, se_type = "stata")
Bacon_Fit1_P5 <- lm_robust(y_p5 ~ Treat + Assembly + Year,
                    data = Bacon_df, cluster = Assembly, se_type = "stata")

Bacon_Fit1_P1などには議会ダミー、年度ダミーが含まれているため、推定結果をひと目で見るのは非常に難しいです。しかし、私たちにとって興味のある推定値は処置変数であるTreatの係数です。{broom}パッケージのtidy()を使うと、係数の推定値と標準誤差、t統計量などを抽出することができ、Treatは2行目に格納されていると考えられます（1行目は切片ですね）。{broom}は{tidyverse}に含まれているため、別途の読み込み作業は不要です。それでは確認してみましょう。

tidy(Bacon_Fit1_P1)[2, ]
##    term   estimate  std.error statistic   p.value    conf.low  conf.high  df
## 2 Treat 0.02153703 0.02784091  0.773575 0.4396403 -0.03319618 0.07627024 399
##   outcome
## 2    y_p1

推定値、標準誤差、t統計量、p値、95%信頼区間などが抽出できました。これは1行のみで構成された1行9列のデータフレームです5。それではここまで推定した5つのモデルから処置効果を抽出し、結合しましょう。データフレームを縦に重ねる際、bind_rows()関数を使います。データフレームの結合については『私たちのR』の第14.5章を参照してください。

Fit1_Coef <- bind_rows(tidy(Bacon_Fit1_P1)[2, ],
                       tidy(Bacon_Fit1_P2)[2, ],
                       tidy(Bacon_Fit1_P3)[2, ],
                       tidy(Bacon_Fit1_P4)[2, ],
                       tidy(Bacon_Fit1_P5)[2, ])

これまでのコードは以下の{purrr}のmap()を使うと、以下のように簡略化することができます。

tibble(Outcome  = paste0("y_p", 1:5),
       Formula1 = paste0(Outcome, "~ Treat + Assembly + Year")) %>%
  mutate(Formula2 = map(Formula1, as.formula),
         Fit      = map(Formula2, ~lm_robust(.x, data = Bacon_df, 
                                             cluster = Assembly, 
                                             se_type = "stata")),
         Tidy     = map(Fit, broom::tidy)) %>%
  select(-(Formula1:Fit)) %>%
  unnest(Tidy) %>%
  filter(term == "Treat")
# 中身の確認
Fit1_Coef
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
term estimate std.error statistic p.value conf.low conf.high df outcome
Treat 0.0215370 0.0278409 0.773575 0.4396403 -0.0331962 0.0762702 399 y_p1
Treat 0.0714665 0.0273206 2.615849 0.0092391 0.0177558 0.1251772 398 y_p2
Treat 0.1186472 0.0278130 4.265884 0.0000249 0.0639671 0.1733273 395 y_p3
Treat 0.1060501 0.0256414 4.135896 0.0000434 0.0556354 0.1564648 385 y_p4
Treat 0.0676668 0.0221121 3.060169 0.0023684 0.0241901 0.1111435 382 y_p5

続いて、それぞれのモデルを識別するための列と統計的有意性を示す列を追加しましょう。また、termやdf、outcome行は不要ですので除外します（あっても無害ですが、画面に出力する際に横長になりますね）。

Fit1_Coef <- Fit1_Coef %>%
    mutate(Model = c("t+1", "t+2", "t+3", "t+4", "t+5"),
           Model = fct_inorder(Model),
           Sig = if_else(conf.low * conf.high > 0, 
                         "Significant", "Insignificant")) %>%
  select(-c(term, df, outcome))
# 中身の確認
Fit1_Coef
estimate std.error statistic p.value conf.low conf.high Model Sig
0.0215370 0.0278409 0.773575 0.4396403 -0.0331962 0.0762702 t+1 Insignificant
0.0714665 0.0273206 2.615849 0.0092391 0.0177558 0.1251772 t+2 Significant
0.1186472 0.0278130 4.265884 0.0000249 0.0639671 0.1733273 t+3 Significant
0.1060501 0.0256414 4.135896 0.0000434 0.0556354 0.1564648 t+4 Significant
0.0676668 0.0221121 3.060169 0.0023684 0.0241901 0.1111435 t+5 Significant

作図の時間です！作図においてある変数の値に応じて色分けを行う場合、幾何オブジェクト（以下の例だと、geom_pointrange()）のaes()内にcolorを指定しますね。色は{ggplot2}が自動的に指定してくれますが、これを自分で調整する場合はscale_color_*()関数群を使用します。詳しい使い方に付いては『私たちのR』の第19.3章を参照してください。

Fit1_Coef %>% 
  ggplot() +
  # y = 0の水平線 (破線)
  geom_hline(yintercept = 0, linetype = 2) +
  # 各モデルの点推定値と信頼区間を出力
  # Sigの値によって色分け
  geom_pointrange(aes(x = Model, y = estimate, 
                      ymin = conf.low, ymax = conf.high,
                      color = Sig), 
                  size = 1.5) +
  # 軸ラベルを設定
  labs(x = "結果変数", y = "DID推定量 w/ 95%信頼区間", color = "") +
  # 色分けを手動で調整
  scale_color_manual(values = c("Significant"   = "black",
                                "Insignificant" = "gray70")) +
  # 日本語フォントの指定
  theme_minimal(base_family = "HiraKakuProN-W3") +
  # 文字サイズ、凡例位置調整
  theme(text = element_text(size = 16),
        legend.position = "bottom")

線形トレンドを考慮したDID推定

これまでのやり方と同じです。ただし、今回はトレンドを表す変数としてTrendが投入されます。このTrendは単体で投入されるのではなく、ユニットのダミー変数であるAssemblyとの交差項で投入します。また、本稿はトレンドが2次曲線であると仮定し、トレンド変数の二乗項を投入しています。ここでもI(Trend^2)で対応しましょう6。

2つの変数を*で繋いだらそれぞれの変数と交差項が同時に投入されますが、:で繋ぐと交差校のみが投入されます。

# Model 2: Diff-in-Diff Estimation with linear trends
Bacon_Fit2_P1 <- lm_robust(y_p1  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_P2 <- lm_robust(y_p2  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_P3 <- lm_robust(y_p3  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_P4 <- lm_robust(y_p4  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_P5 <- lm_robust(y_p5  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")

以下の作業は線形トレンドを仮定しなかったモデルの可視化と同じ手順となります。

Fit2_Coef <- bind_rows(tidy(Bacon_Fit2_P1)[2, ],
                       tidy(Bacon_Fit2_P2)[2, ],
                       tidy(Bacon_Fit2_P3)[2, ],
                       tidy(Bacon_Fit2_P4)[2, ],
                       tidy(Bacon_Fit2_P5)[2, ])
Fit2_Coef <- Fit2_Coef %>%
    mutate(Model = c("t+1", "t+2", "t+3", "t+4", "t+5"),
           Model = fct_inorder(Model),
           Sig = if_else(conf.low * conf.high > 0, 
                         "Significant", "Insignificant")) %>%
  select(-c(term, df, outcome))
# 中身の確認
Fit2_Coef
estimate std.error statistic p.value conf.low conf.high Model Sig
0.0488728 0.0232529 2.101793 0.0361978 0.0031593 0.0945862 t+1 Significant
0.0533259 0.0241351 2.209470 0.0277107 0.0058776 0.1007742 t+2 Significant
0.0616731 0.0254269 2.425502 0.0157345 0.0116840 0.1116622 t+3 Significant
0.0521781 0.0249654 2.090019 0.0372717 0.0030926 0.1012637 t+4 Significant
0.0261356 0.0237592 1.100020 0.2720160 -0.0205796 0.0728508 t+5 Insignificant 
Fit2_Coef %>% 
  ggplot() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_pointrange(aes(x = Model, y = estimate, 
                      ymin = conf.low, ymax = conf.high,
                      color = Sig), 
                  size = 1.5) +
  labs(x = "結果変数", y = "DID推定量 w/ 95%信頼区間", color = "") +
  scale_color_manual(values = c("Significant"   = "black",
                                "Insignificant" = "gray70")) +
  theme_minimal(base_family = "HiraKakuProN-W3") +
  theme(text = element_text(size = 16),
        legend.position = "bottom")

頑健性の確認

頑健性のチェックは様々な方法が考えられますが、今回は本稿の著者たちが行った方法で確認してみます。それは\(t\)年や\(t-1\)なども結果変数として用いることです。もし、処置によって過去の補助金額が増額されたら、\(t+1\)期以降に見られる補助金増額は処置による効果とは言いにくくなりますね。

\(t\)年の補助金額はy_0、\(t-1\)年のそれはy_m1です。今回はトレンド変数 (の二乗項)ありのモデルを使います。

Bacon_Fit2_M5 <- lm_robust(y_m5  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_M4 <- lm_robust(y_m4  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_M3 <- lm_robust(y_m3  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_M2 <- lm_robust(y_m2  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_M1 <- lm_robust(y_m1  ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Bacon_Fit2_0  <- lm_robust(y_0   ~ Treat + Assembly + Year +
                             Assembly:I(Trend^2), data = Bacon_df,
                           cluster = Assembly, se_type = "stata")
Fit2_Coef2 <- bind_rows(tidy(Bacon_Fit2_M5)[2, ],
                        tidy(Bacon_Fit2_M4)[2, ],
                        tidy(Bacon_Fit2_M3)[2, ],
                        tidy(Bacon_Fit2_M2)[2, ],
                        tidy(Bacon_Fit2_M1)[2, ],
                        tidy(Bacon_Fit2_0)[2, ])
Fit2_Coef2 <- Fit2_Coef2 %>%
  mutate(Model = c("t-5", "t-4", "t-3", "t-2", "t-1", "t"),
         Model = fct_inorder(Model),
         Sig = if_else(conf.low * conf.high > 0, 
                       "Significant", "Insignificant")) %>%
  select(-c(term, df, outcome))
# 中身の確認
Fit2_Coef2
estimate std.error statistic p.value conf.low conf.high Model Sig
0.0224672 0.0235360 0.9545907 0.3403821 -0.0238075 0.0687420 t-5 Insignificant
-0.0285913 0.0212211 -1.3473037 0.1786687 -0.0703141 0.0131315 t-4 Insignificant
-0.0028826 0.0216179 -0.1333416 0.8939904 -0.0453819 0.0396167 t-3 Insignificant
0.0085610 0.0185290 0.4620323 0.6443100 -0.0278657 0.0449877 t-2 Insignificant
0.0311195 0.0183498 1.6958998 0.0906849 -0.0049550 0.0671940 t-1 Insignificant
0.0393396 0.0234982 1.6741546 0.0948839 -0.0068561 0.0855352 t Insignificant

続いて、Fit2_Coefと結合します。また、できれば可視化の際、横軸を\(t-5\)から\(t+5\)の順で並べたいので、もう一回、fct_inorder()を使って順番を指定しておきます。

Robust_Coef <- bind_rows(Fit2_Coef2, Fit2_Coef)
# 中身をチェック
Robust_Coef
estimate std.error statistic p.value conf.low conf.high Model Sig
0.0224672 0.0235360 0.9545907 0.3403821 -0.0238075 0.0687420 t-5 Insignificant
-0.0285913 0.0212211 -1.3473037 0.1786687 -0.0703141 0.0131315 t-4 Insignificant
-0.0028826 0.0216179 -0.1333416 0.8939904 -0.0453819 0.0396167 t-3 Insignificant
0.0085610 0.0185290 0.4620323 0.6443100 -0.0278657 0.0449877 t-2 Insignificant
0.0311195 0.0183498 1.6958998 0.0906849 -0.0049550 0.0671940 t-1 Insignificant
0.0393396 0.0234982 1.6741546 0.0948839 -0.0068561 0.0855352 t Insignificant
0.0488728 0.0232529 2.1017929 0.0361978 0.0031593 0.0945862 t+1 Significant
0.0533259 0.0241351 2.2094697 0.0277107 0.0058776 0.1007742 t+2 Significant
0.0616731 0.0254269 2.4255022 0.0157345 0.0116840 0.1116622 t+3 Significant
0.0521781 0.0249654 2.0900194 0.0372717 0.0030926 0.1012637 t+4 Significant
0.0261356 0.0237592 1.1000203 0.2720160 -0.0205796 0.0728508 t+5 Insignificant 
Robust_Coef %>%
  ggplot() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_pointrange(aes(x = Model, y = estimate, 
                      ymin = conf.low, ymax = conf.high,
                      color = Sig), 
                  size = 1) +
  labs(x = "結果変数", y = "DID推定量 w/ 95%信頼区間", color = "") +
  scale_color_manual(values = c("Significant"   = "black",
                                "Insignificant" = "gray70")) +
  theme_minimal(base_family = "HiraKakuProN-W3") +
  theme(text = element_text(size = 16),
        legend.position = "bottom")

Fouirnaies and Mutlu-Eren (2015)のFigure 3に比べると結果が若干異なります。なぜなら、Fouirnaies and Mutlu-Eren (2015)で推定したモデルとここで行った分析のモデルが一致しないからです。ここでは投入しなかった人口や選挙サイクルのダミー変数がFouirnaies and Mutlu-Eren (2015)では投入されています。また、ここで提供しているデータはFouirnaies and Mutlu-Eren (2015)と同じやり方でクリーニングされていない点も一つの原因です。それでもFouirnaies and Mutlu-Eren (2015)のFigure 3の非常に似たような結果が得られました。興味のある方は著者たちが提供している再現用データセットをダウンロードし、取り組んでみましょう。

Synthetic Control Method

{gsynth}パッケージの使い方

それでは、差分の差分法の応用としてSynthetic Control Method (SCM) とCausalImpactについて解説します。ここではパッケージの使い方に重点を起きながら説明します。

まずは、SCMからです。Synthetic Control MethodをRに実装したパッケージは{Synth}、{gsynth}、{bpCausal}などがあります。ここでは使い方が最も簡単な{gsynth}パッケージについて説明します。また、後半では{Synth}を利用したAbadie et al. (2015)の再現コードも掲載します。

ここで使用するのは歴代参院選（第7回以降）における都道府県の投票率です7。処置変数は第24回参院選から導入された「合区」です。具体的には鳥取選挙区と島根選挙区が「鳥取・島根選挙区」に、徳島選挙区と高知選挙区が「徳島・高知選挙区」として合区されたか否かです。合区によって自分の1票の価値がほぼ半分に下がり、投票参加を妨げたのではないかという議論もありましたが、本当でしょうか。

{gsynth}と{panelView}パッケージを読み込み、歴代参院選の都道府県別投票率データ（Day3_Data6.csv）を読み込みます。パッケージ読み込みの際、{pacman}のp_load()を使うと、インストールされていないパッケージをインストールしてから読み込んでくれるので、大変便利です。

pacman::p_load(gsynth, panelView)

upper_df <- read_csv("Data/Day3_Data6.csv")
## Rows: 891 Columns: 9
## ─ Column specification ────────────────────────────
## Delimiter: ","
## chr (2): PrefName_E, PrefName_J
## dbl (7): Election, PrefCode, Eligible, Voters, EffVote, Magnitude, Candidates
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

記述統計を確認してみましょう。

view(dfSummary(upper_df))
No Variable Stats / Values Freqs (% of Valid) Graph Valid Missing
1 Election [numeric] Mean (sd) : 16 (5.5) min < med < max: 7 < 16 < 25 IQR (CV) : 10 (0.3) 19 distinct values 891 (100.0%) 0 (0.0%)
2 PrefCode [numeric] Mean (sd) : 23.9 (13.5) min < med < max: 1 < 24 < 47 IQR (CV) : 24 (0.6) 47 distinct values 891 (100.0%) 0 (0.0%)
3 PrefName_E [character] 1. Aichi 2. Akita 3. Aomori 4. Chiba 5. Ehime 6. Fukui 7. Fukuoka 8. Fukushima 9. Gifu 10. Gunma [ 37 others ]
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
701(78.7%)
 891 (100.0%) 0 (0.0%)
4 PrefName_J [character] 1. 愛知県 2. 愛媛県 3. 茨城県 4. 岡山県 5. 岩手県 6. 岐阜県 7. 宮崎県 8. 宮城県 9. 京都府 10. 熊本県 [ 37 others ]
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
19(2.1%)
701(78.7%)
 891 (100.0%) 0 (0.0%)
5 Eligible [numeric] Mean (sd) : 1917557 (1842564) min < med < max: 364184 < 1189566 < 11396789 IQR (CV) : 1108658 (1) 891 distinct values 891 (100.0%) 0 (0.0%)
6 Voters [numeric] Mean (sd) : 1138294 (1016940) min < med < max: 237073 < 764645 < 6415327 IQR (CV) : 639952 (0.9) 889 distinct values 891 (100.0%) 0 (0.0%)
7 EffVote [numeric] Mean (sd) : 1095111 (985815.3) min < med < max: 228792 < 729981 < 6223547 IQR (CV) : 610682 (0.9) 891 distinct values 891 (100.0%) 0 (0.0%)
8 Magnitude [numeric] Mean (sd) : 1.6 (0.9) min < med < max: 1 < 1 < 6 IQR (CV) : 1 (0.5)
1:517(58.0%)
2:258(29.0%)
3:78(8.8%)
4:33(3.7%)
5:3(0.3%)
6:2(0.2%)
 891 (100.0%) 0 (0.0%)
9 Candidates [numeric] Mean (sd) : 5.4 (5.4) min < med < max: 2 < 4 < 72 IQR (CV) : 2 (1) 30 distinct values 891 (100.0%) 0 (0.0%)

Generated by summarytools 0.9.9 (R version 4.0.4)
2021-08-04

変数名 説明 備考
Election 選挙 第XX回参議院議員通常選挙
PrefCode 都道府県 JIS規格コード
PrefName_E 都道府県 英語
PrefName_J 都道府県 日本語
Eligible 有権者数
Voters 投票者数
EffVote 有効投票数
Magnitude 定数
Candidates 候補者数

合区された選挙区における定数および候補者数の扱いですが、たとえば島根選挙区の場合、合区前の定数は1です。合区後の鳥取・島根選挙区も定数1です。この場合、合区後の島根「県」の定数を1にするか0.5にするか、そして候補者数も合区における候補者数にすべきか、2に割るかが問題となりますが、ここでは定数を1とし、候補者数も合区における候補者数にしました。このようなコーディングが実証上、正しいかどうかは分かりませんが、実習の場面では問題ないでしょう。

それでは、次はデータのクリーニングします。ここでは4つの作業を行います。

  1. Treatという変数を作成する。PrefNameが「鳥取県」、「島根県」、「徳島県」、「高知県」であり、Electionが24以上なら1、それ以外のケースは0とする。
  2. Turnoutという変数を作成する。投票者数（Voters）を有権者数（Eligible）で割り、100を掛ける。
  3. Spoiltという変数を作成する。有効投票数（EffVote）を有権者数（Eligible）で割った値を1から引き、100を掛ける。
  4. 都道府県名（PrefName_E）をfactor変数に変換し、PrefName_Fという列として追加します。水準（levels）の順番はデータの出現順とする。
upper_df <- upper_df %>%
  mutate(Treat      = if_else(PrefName_E %in% c("Tottori", "Shimane",
                                                "Tokushima","Kochi") &
                                Election >= 24, 1, 0),
         Turnout    = Voters / Eligible * 100,
         Spoilt     = (1 - (EffVote / Voters)) * 100,
         PrefName_F = fct_inorder(PrefName_E))

本格的な分析に入る前に、各ケースの処置有無を可視化してみましょう。先ほど読み込みました{panelView}パッケージを使います。第一引数は応答変数 ~ 処置変数であり、data引数にデータフレームのオブジェクト名を当てます。また、indexにはユニットと時間を表す変数名を長さ2のcharacter型ベクトルとして指定します。ここでは"PrefName_E"と"Election"ですね。最後に、pre.post = TRUEを指定すると、処置前後を色分けしてくれるので、見やすくなります。

panelView(Turnout ~ Treat, data = upper_df, 
          index = c("PrefName_F", "Election"), pre.post = TRUE,
          xlab = "Election", ylab = "Prefecture")

これを見ると処置を受けるケースが4県であり、どれも第24回参院選から処置を受けることが分かります。また、沖縄県の場合、第7・8回参院選のデータが欠損していることが分かります。

次は投票率の変化を時系列的に示してみましょう。使い方は先ほどとほぼ同じですが、panelView()内にtype = "outcome"引数を追加します。

panelView(Turnout ~ Treat, data = upper_df, 
          index = c("PrefName_F", "Election"), type = "outcome")

これを見ると、都道府県内における投票率の変化にはバラツキがありますが、傾向としてはかなり似通っていることが分かります。また、処置後の変化（赤い線）を見ると、他の都道府県よりかなり落ちているようにも見えますね。これは合区によって何かの変化が生じた可能性があることを示唆しています。

それではSCMを実行してみましょう。基本的な使い方はlm()関数に近いですが、それでも必要な引数がそれなりにあります。ここではその一部を紹介します。

  • 第1引数は応答変数 ~ 処置変数 + 統制変数1 + 統制変数2 + ...であり、統制変数は必須ではありません。
  • dataにはデータフレームのオブジェクト名を指定します。
  • indexはpanelView()と同じですが、一つ注意が必要です。それは、ユニットを表す変数をこれまでは"PrefName_F"にしましたが、ここでは"PrefName_E"にすることです。まだ開発途上のパッケージということもあり、ユニット変数がfactor型だと、正しく推定できません。一方、PrefName_Eはfactor型でなくcharacter型ですので、問題ございません。また、numeric型であるPrefCodeを指定してもOKです。
  • forceは固定効果をユニットレベルにするか、時間レベルにするか、両方にするか、あるいはなしにするかです。既定値はユニットレベル（"unit"）ですが、ここでは両方（"two-way"）にしてみます。
  • seは標準誤差を計算するか否かであり、既定値はFALSEです。ここでは標準誤差も計算するためにTRUEにしました。標準誤差を計算しない場合、計算が早く終わります。
  • inferenceは推定方法ですが、既定値はノンパラメトリック推定（"nonparametric"）です。ただし、処置ケースが少ない場合、パラメトリック推定が推奨されているため、ここでは"parametric"を指定します。
  • nbootsは標準誤差を計算する際に使用されるブートストラッピングの回数です。既定値は200ですが、ここでは1000回にします。
  • parallelは並列計算の有無ですが、既定値はTRUEで、そのままで問題ありません。ただし、R Markdownに入れる場合、並列計算ができないため、ここではFALSEにしておきます。
synth_fit <- gsynth(Turnout ~ Treat + Magnitude + Candidates, 
                    data = upper_df, index = c("PrefName_E", "Election"), 
                    force = "two-way", se = TRUE, inference = "parametric",
                    nboots = 1000, parallel = FALSE)
## Cross-validating ... 
##  r = 0; sigma2 = 8.35872; IC = 2.13975; MSPE = 12.90393
##  r = 1; sigma2 = 6.09626; IC = 2.32584; MSPE = 11.26215
##  r = 2; sigma2 = 5.09782; IC = 2.63224; MSPE = 9.68534
##  r = 3; sigma2 = 4.28768; IC = 2.92798; MSPE = 9.22076*
##  r = 4; sigma2 = 3.86785; IC = 3.27729; MSPE = 9.40547
##  r = 5; sigma2 = 3.44595; IC = 3.59771; MSPE = 11.04103
## 
##  r* = 3
## 
## 
Simulating errors .............
Bootstrapping ...
## ..........
print(synth_fit)
## Call:
## gsynth.formula(formula = Turnout ~ Treat + Magnitude + Candidates, 
##     data = upper_df, index = c("PrefName_E", "Election"), force = "two-way", 
##     se = TRUE, nboots = 1000, inference = "parametric", parallel = FALSE)
## 
## Average Treatment Effect on the Treated:
##  ATT.avg  S.E. CI.lower CI.upper p.value
##   -6.755 1.505   -9.499   -3.475       0
## 
##    ~ by Period (including Pre-treatment Periods):
##         ATT   S.E. CI.lower CI.upper p.value n.Treated
## 1   0.06989 1.1453  -1.7736  2.62214   0.708         0
## 2  -0.68512 1.1783  -2.6418  1.92466   0.756         0
## 3   3.30567 1.0448   1.5806  5.58852   0.000         0
## 4   1.37507 1.0935  -0.8522  3.41191   0.220         0
## 5  -1.42366 1.1102  -3.7058  0.58751   0.132         0
## 6  -1.14623 0.9409  -3.4557  0.16800   0.102         0
## 7  -2.62647 1.1161  -4.6741 -0.29892   0.032         0
## 8  -1.60347 1.0930  -4.1767  0.15326   0.082         0
## 9   1.22010 1.0217  -0.8904  3.21370   0.262         0
## 10 -0.86152 1.0392  -2.6675  1.51308   0.580         0
## 11  3.51810 1.5325   0.9092  6.99112   0.014         0
## 12 -0.94563 1.0699  -2.9268  1.33867   0.456         0
## 13  0.34662 0.9317  -1.5177  2.08426   0.732         0
## 14 -0.37951 0.9785  -2.4387  1.37640   0.628         0
## 15  0.26459 0.9162  -1.6980  1.98752   0.844         0
## 16  1.60321 0.9221  -0.4818  3.07703   0.174         0
## 17 -2.03164 0.9435  -3.8030 -0.05776   0.044         0
## 18 -6.58463 1.3882  -9.3547 -3.70984   0.000         4
## 19 -6.92573 1.9747 -10.4256 -2.55373   0.004         4
## 
## Coefficients for the Covariates:
##                beta    S.E. CI.lower CI.upper p.value
## Magnitude   2.31026 0.77861   0.9812  4.04461   0.000
## Candidates -0.07379 0.04336  -0.1513  0.02325   0.126

平均値な処置効果（ATT）は-6.755であり、統計的にも有意な結果が得られました。これは合区が行われた選挙区は、もし合区しなかった場合の投票率に比べ、約6.755%ポイント低いことを意味します。また、~ by Period以下の欄では各選挙ごとのATTが表示されます。18と19が合区以降の時期であり、それぞれ約-6.585%p、-6.926%pです。上記の-6.755%pはこの2つの数値の平均値です。ちなみに処置群の県が4つにも関わらず、ATTが一つだけ表示されるのは、4つの平均値を出しているからです。それぞれの県ごとに効果を確認する方法は後ほど解説します。とりあえず、この結果から合区は投票率を下げたという解釈ができるでしょう。

これを可視化するためにはplot()関数を使います。type = "gap"を指定すると、もし処置を受けなかった場合の4県の平均値と実際の平均値が折れ線グラフで出力されます。

plot(synth_fit, type = "gap", main = "Estimated ATT", 
     xlab = "Election (0 = 1990 Election)", ylab = "Gap (%p)")

処置を受ける前は架空の4県とその他の43都道府県にはあまり差がありませんでしたが、処置を受けてから差が広がることが分かります。差分でなく、架空の4県と実際の4県のトレンドを見るためにはtype = "counterfactual"、もしくはtype = "ct"を指定します。

plot(synth_fit, type = "counterfactual", 
     xlab = "Election", ylab = "Turnout (%)")

処置を受けたのは4県ですが、平均値を取り、線が1本に圧縮しました。ここからも処置前は架空の4県と実際の4県はほぼ同じトレンドを示していますが、処置後は傾向が変わりましたね。

もし、4つの都道府県を個別に示したい場合はどうすれば良いでしょうか。{gsynth}から得られたオブジェクトには架空の鳥取、島根、徳島、高知のデータが個別に格納されているため、それを抽出して使います。ここではコードのみをお見せします。

# 処置群の実際のデータは「$Y.tr」で抽出可能。matrix型として格納されている。
Treat_df <- synth_fit$Y.tr %>%
  as_tibble() %>%
  mutate(Election = 7:25, .before = Kochi) %>%
  pivot_longer(cols      = Kochi:Tottori,
               names_to  = "Pref",
               values_to = "Turnout")

# 処置群が処置を受けなかった場合の架空ののデータは「$Y.ct」で抽出可能。
# matrix型として格納されている。
Counter_df <- synth_fit$Y.ct %>%
  as_tibble() %>%
  mutate(Election = 7:25, .before = Kochi) %>%
  pivot_longer(cols      = Kochi:Tottori,
               names_to  = "Pref",
               values_to = "Turnout")

bind_rows(list("Actual" = Treat_df, "Counterfactual" = Counter_df),
          .id = "Type") %>%
  mutate(Type = factor(Type, levels = c("Counterfactual", "Actual"))) %>%
  ggplot() +
  geom_vline(xintercept = 23, linetype = 2) +
  geom_line(aes(x = Election, y = Turnout, linetype = Type, color = Type),
            size = 1) +
  scale_color_manual(values = c("Actual" = "black", 
                                "Counterfactual" = "orange")) +
  scale_linetype_manual(values = c("Actual" = 1, 
                                   "Counterfactual" = 2)) +
  guides(linetype = "none") +
  labs(x = "Election", y = "Turnout (%)", color = "") +
  facet_wrap(~Pref, ncol = 2) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

# Gapのデータは「$est.ind」から抽出可能。3次元array型として格納されている。
bind_rows(list("Tottori"   = as_tibble(synth_fit$est.ind[, , "Tottori"]),
               "Shimane"   = as_tibble(synth_fit$est.ind[, , "Shimane"]),
               "Tokushima" = as_tibble(synth_fit$est.ind[, , "Tokushima"]),
               "Kochi"     = as_tibble(synth_fit$est.ind[, , "Kochi"])),
          .id = "Pref") %>%
  mutate(Election = rep(7:25, 4), .before = Pref) %>%
  mutate(Pref = fct_inorder(Pref)) %>%
  ggplot() +
  geom_hline(yintercept = 0, linetype = 2) +
  geom_vline(xintercept = 23, linetype = 2) +
  geom_ribbon(aes(x = Election, ymin = CI.lower, ymax = CI.upper),
              alpha = 0.25) +
  geom_line(aes(x = Election, y = EFF)) +
  labs(x = "Election", y = "Gap (%p)") +
  facet_wrap(~Pref, ncol = 2) +
  theme_bw(base_size = 12) +
  theme(legend.position = "bottom")

個別に見ると、投票率を下げる効果が会ったのは鳥取県と徳島県で、島根県では効果が確認されていません。また、高知の場合、合区直後のみにおいて効果が観察されます。

SCMは統制群から合成された架空の鳥取、島根、徳島、高知を作成し、こちらを実際の統制群として用います。その際に43都道府県を重み付け合成するわけですが、それぞれの重みはいくつでしょうか。それを確認するためにはgsynth()から得られたオブジェクトからwgt.impliedを抽出します。

synth_fit$wgt.implied
##                 Kochi     Shimane  Tokushima     Tottori
## Aichi     -0.47478433 -3.11286790 -3.0401794  1.46035320
## Akita     -0.23531295 -0.52943972 -0.9885995  0.44729458
## Aomori    -0.33133036 -1.42351775 -1.5004957  0.63697503
## Chiba     -0.10072384 -2.72848780 -1.9602985  1.06654903
## Ehime      0.36693421  0.06277326  0.4126363  0.06224253
## Fukui      0.22372066  6.13863419  5.3608227 -3.11199739
## Fukuoka   -0.35968294  0.85837901  0.4556590 -0.60296813
## Fukushima  0.92413466 -0.10315708  0.7606631  0.33627698
## Gifu       0.34763407  0.40058919  0.6455275 -0.09395036
## Gunma      0.81795260  2.77601904  3.0196358 -1.14178187
## Hiroshima  0.12536302  2.10005857  2.1181854 -1.22844930
## Hokkaido  -1.24999673 -4.38247164 -5.0753485  2.07357037
## Hyogo     -0.53000971 -3.93684511 -3.8384699  1.91187976
## Ibaraki    0.39650251 -0.57436139  0.7862381 -0.25676441
## Ishikawa  -0.10753601  0.61996133  0.3688639 -0.31437700
## Iwate     -0.61873142 -3.91679568 -4.3203523  2.18592914
## Kagawa     0.30182287  4.75543151  4.3745006 -2.44893198
## Kagoshima  0.15709875 -0.11220145 -0.2996091  0.36706375
## Kanagawa  -0.25513789 -4.69770311 -4.0302773  2.24448298
## Kumamoto  -0.42759500  2.16670497  1.1126506 -1.01823375
## Kyoto     -0.55809809 -5.22670116 -5.2072717  2.76546457
## Mie       -0.23120385 -1.91335402 -1.7577212  0.87484030
## Miyagi     0.67548370  0.16388736  1.1355756 -0.16082668
## Miyazaki   0.17920426  7.55529258  6.6907730 -3.99468602
## Nagano     0.18001918 -1.38692919 -0.9438182  0.73188001
## Nagasaki  -0.50914219  2.09290651  1.3335133 -1.27192928
## Nara      -0.15351003 -2.78581544 -2.3565366  1.30382970
## Nigata     0.49147410 -2.78204340 -1.4733561  1.28378025
## Oita      -0.17349464  2.79744500  1.7438750 -1.16909040
## Okayama    0.07652751  0.96627526  1.3455797 -0.82350635
## Okinawa    1.01863146  0.42161368  0.6455559  0.57295931
## Osaka     -0.44577523 -4.47231131 -4.6474027  2.54090911
## Saga       0.12277513  6.48755504  5.9408655 -3.60701753
## Saitama    0.05216851 -1.75601014 -0.8541794  0.49812787
## Shiga     -0.11386770  0.33370104 -0.2356226  0.09945418
## Shizuoka   0.91971452  1.54285585  2.7581497 -0.98006477
## Tochigi    0.52297922  1.00450166  1.7008602 -0.64235400
## Tokyo     -0.38924399 -5.94743192 -5.6164115  3.15967284
## Toyama     0.16744143  6.06393887  5.6432566 -3.38672772
## Wakayama  -0.31260079 -1.62066751 -1.8282660  0.87875300
## Yamagata   0.39664997 -0.94612490 -0.5127014  0.67261421
## Yamaguchi -0.75068452  1.76262527  0.1130550 -0.63602344
## Yamanashi -0.13577011  3.28408845  2.0204754 -1.28522233

43行4列のデータが得られました。たとえば、架空の高知県を作成する時に、北海道と沖縄県もデータが多く参照されていることが分かります。

Diff-in-Diffで推定してみたら?

upper_df %>%
  # トレンド変数を作成する
  # 実はElection変数は1単位で増加しているため、作らなくても良い
  mutate(Trend = Election - 6) %>%
  lm_robust(Turnout ~ Treat + Magnitude + Candidates + 
              PrefName_F + factor(Election) + PrefName_F:Trend, 
            data = ., cluster = PrefName_F, se_type = "stata") %>%
  broom::tidy() %>%
  filter(term == "Treat")
## 1 coefficient  not defined because the design matrix is rank deficient
##    term  estimate std.error statistic      p.value  conf.low conf.high df
## 1 Treat -5.317369 0.8201135 -6.483699 5.435995e-08 -6.968172 -3.666566 46
##   outcome
## 1 Turnout

合区が投票率に与えた影響は約-5.317%p [-6.968, -3.667]と推定されました。一般化SCMの推定値より約1.4%p小さいですね。

{Synth}を用いたAbadie et al. (2015)の再現

RでSCMを実装する場合、広く使われてきたパッケージとして{Synth}があります。こちらは使い方が非常に難しいですが、ここではAbadie et al. (2015)の結果を再現してみましょう8。まずは、{Synth}パッケージと実習用データを読み込みましょう。こちらはAbadie et al. (2015)9のリプリケーションデータをcsvに変換したものとなります。

library(Synth)
## ##
## ## Synth Package: Implements Synthetic Control Methods.
## ## See http://www.mit.edu/~jhainm/software.htm for additional information.
# read_csvでなく、read.csvを使用
# R 4.0.0未満の場合、stringsAsFactors = FALSEを付けること
Synth_Data <- read.csv("Data/Day3_Data3.csv")

データの中身を確認してみましょう。

library(DT)
datatable(Synth_Data)

データの詳細はAbadie et al. (2015)の本文中付録を参照してください。ここで必要なものは以下の変数です。

変数 説明
index 国ID
country 国名
year 年
gdp 一人当たりGDP
infrate インフレーション率
trade 貿易依存度
schooling 25歳以上で中等教育以上修了者の割合
invest70 国内総投資の割合
industry 高付加価値産業のシェア

まずは、可視化をしてみましょう。横軸は年、縦軸は一人当たりGDPです。西ドイツ以外の国は平均を計算し、一本の線としてまとめます。また、ドイツが統一した1990年に垂直線を引きます。

#quartz(type = "pdf", file = "../Figs/SCM2.pdf", width = 8, height = 4)
Synth_Data %>%
  mutate(WG = if_else(country == "West Germany", 
                      "西ドイツ", "西ドイツ以外のOECD加盟国")) %>%
  group_by(WG, year) %>%
  summarise(GDP     = mean(gdp),
            .groups = "drop") %>%
  ggplot() +
  geom_vline(xintercept = 1990, linetype = 2) +
  geom_line(aes(x = year, y = GDP, color = WG), size = 1) +
  annotate("text", x = 1990, y = 24000, label = "ドイツ統一 → ",
           hjust = 1, family = "HiraKakuProN-W3", size = 5) +
  labs(x = "年", y = "GDP (米ドル)", 
       color = "", linetype = "") +
  theme_minimal(base_family = "HiraKakuProN-W3") +
  theme(legend.position = "bottom",
        text = element_text(size = 16))

#dev.off()

赤い線が西ドイツ、青がその他のOECD加盟国の平均GDPです。全体的な傾向は同じですが、80年代以降はトレンドが並行するというよりも、どんどん拡大しているように見えます。

それでは、SCMの実装について説明します。{Synth}パッケージのメインは、データの準備に使うdataprep()とSCMの推定に使うsynth()です。snyth()は比較的に使いやすいですが、問題はdataprep()です。ここでは、個々のパラメーターについて詳細に解説するより、コードの解説を中心に進めます。

Train_Data <- dataprep(
  # データを指定
  foo           = Synth_Data,
  # トレンドを予測する際に使用する共変量を指定
  # ここではgdp, infrate, tradeを使います
  predictors    = c("gdp", "infrate", "trade"),
  # 結果変数を指定
  dependent     = "gdp",
  # 観察対象を指定
  # ここではindexとcountry列がありますが、ここはnumeric型である必要があり、
  # indexを指定します。countryは後ほど使います
  unit.variable = "index",
  # 時間を表す変数を指定
  time.variable = "year",
  # ここではやや特殊な共変量を指定します
  # たとえば、1971年から1980年までのindustryの平均値などです
  special.predictors = list(
    list("industry",  1971:1980,    c("mean")),
    list("schooling", c(1970,1975), c("mean")),
    list("invest70",  1980,         c("mean"))
  ),
  # 処置変数を受ける観察対象の番号を指定します
  # West Germanyのindex番号は7です
  treatment.identifier  = 7,
  # 処置群以外の観察対象の番号を指定します
  controls.identifier   = unique(Synth_Data$index)[-7],
  # 処置前の期間を指定します
  # year変数が1971から1980までの場合を処置前とします（後で説明）
  time.predictors.prior = 1971:1980,
  # 損失関数を最小化する範囲を指定します
  # 基本的には処置前の任意の時点〜処置を受ける時点に設定します
  # ここでは処置前の10年分とします
  time.optimize.ssr     = 1981:1990,
  # 国名が格納されている列名を指定します
  unit.names.variable   = "country",
  # 可視化の際、横軸の範囲を指定します
  time.plot             = 1960:2003
  )

非常に使いにくい関数ですが、それでも自分で直接SCMを実装するよりは遥かに楽です。SCMを行う前に、いくつか説明をしておきます。

まず、time.predictors.prior引数ですが、ここでは1971:1980に指定しました。「ドイツ統一は1990年だから1971:1990にすべきでは？」と思うのが普通でしょう。それが正しいです。ただし、synth()関数を使う場合はいくつかのパラメーターを指定する必要があり、それはsynth()関数が自動的に指定してくれます。しかし、データを100%使ってそのパラメーターを決める場合、過適合（overfitting）が生じる可能性があります。まずは、架空の西ドイツを作るのが目的であるため、ここでは全期間でなく、処置前の一部の期間のみで架空の西ドイツを作成し、そこから得られたパラメーターをもう一回使います。だから、ここではまず1971年から1980年までのデータを使います。

それではSCMを行います。先ほどdataprep()から作成したオブジェクトをsynth()関数に渡すわけです。

Train_Fit <- synth(Train_Data)
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  searching for synthetic control unit  
##  
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 4556.721 
## 
## solution.v:
##  0.5598212 0.04804541 0.099915 0.004093611 0.07931743 0.2088073 
## 
## solution.w:
##  0.132081 5.52698e-05 0.5028072 3.72175e-05 0.0001515816 4.09008e-05 3.65981e-05 0.00012152 3.73332e-05 0.1693102 0.1488914 2.72698e-05 2.61162e-05 2.71489e-05 0.04555833 0.0007909723

ここまでの結果を見ることもできますが、先が長いのでとりあえず、本分析に入ります。まずは、本分析のためのデータを用意します。

Main_Data <- dataprep(
  foo           = Synth_Data,
  predictors    = c("gdp","trade","infrate"),
  dependent     = "gdp",
  unit.variable = "index",
  time.variable = "year",
  special.predictors = list(
    list("industry",  1981:1990,    c("mean")), # <<
    list("schooling", c(1980,1985), c("mean")), # <<
    list("invest80",  1980,         c("mean"))  # <<
  ),
  treatment.identifier  = 7,
  controls.identifier   = unique(Synth_Data$index)[-7],
  time.predictors.prior = 1981:1990, # <<
  time.optimize.ssr     = 1960:1989, # <<
  unit.names.variable   = "country",
  time.plot             = 1960:2003
  )
## 
##  Missing data: treated unit; special predictor: special.industry.1981.1990 ; for period: 1990 
##  We ignore (na.rm = TRUE) all missing values for predictors.op.

SCMに使うデータの期間が変わるため、いくつかの行を修正する必要があります。それではSCMを行います。ここではvパラメーターを直接指定します。custom.v = as.numeric(Train_Fit$solution.v)を指定すると、Train_Fitで決められたvパラメーターを本分析のモデルで使うことができます。

Main_Fit <- synth(Main_Data,
                  custom.v = as.numeric(Train_Fit$solution.v))
## 
## X1, X0, Z1, Z0 all come directly from dataprep object.
## 
## 
## **************** 
##  optimization over w weights: computing synthtic control unit 
##  
## 
## 
## **************** 
##  v weights supplied manually: computing synthtic control unit 
##  
## 
## 
## **************** 
## **************** 
## **************** 
## 
## MSPE (LOSS V): 14035.22 
## 
## solution.v:
##  0.5598212 0.04804541 0.099915 0.004093611 0.07931743 0.2088073 
## 
## solution.w:
##  0.2159589 5.2323e-06 0.4017819 1.05206e-05 8.1855e-06 6.9831e-06 4.1299e-06 0.1073192 5.4163e-06 0.1134424 0.1614399 2.2224e-06 2.4139e-06 4.2491e-06 4.7768e-06 3.5972e-06

可視化の前に、バランスチェックと各国の重みを計算してみましょう。ここではsynth.tab()関数を使います。必要な引数はSCMに使用したデータ（dataprep.res =）、SCM推定結果（synth.res =）の2つです。

Synth_Table <- synth.tab(dataprep.res = Main_Data,
                         synth.res    = Main_Fit)

Synth_Table
## $tab.pred
##                               Treated Synthetic Sample Mean
## gdp                         15808.900 15801.364   13669.381
## trade                          56.778    57.346      59.831
## infrate                         2.595     3.434       7.617
## special.industry.1981.1990     34.538    34.385      33.794
## special.schooling.1980.1985    55.500    54.949      38.659
## special.invest80.1980          27.018    27.043      25.895
## 
## $tab.v
##                             v.weights
## gdp                             0.560
## trade                           0.048
## infrate                         0.100
## special.industry.1981.1990      0.004
## special.schooling.1980.1985     0.079
## special.invest80.1980           0.209
## 
## $tab.w
##    w.weights  unit.names unit.numbers
## 1      0.216         USA            1
## 2      0.000          UK            2
## 3      0.402     Austria            3
## 4      0.000     Belgium            4
## 5      0.000     Denmark            5
## 6      0.000      France            6
## 8      0.000       Italy            8
## 9      0.107 Netherlands            9
## 10     0.000      Norway           10
## 12     0.113 Switzerland           12
## 14     0.161       Japan           14
## 16     0.000      Greece           16
## 18     0.000    Portugal           18
## 19     0.000       Spain           19
## 20     0.000   Australia           20
## 21     0.000 New Zealand           21
## 
## $tab.loss
##           Loss W   Loss V
## [1,] 0.003111837 14035.22

論文で主に使うのは$tab.predと$tab.wです。

続いてSCMでよく見る図を可視化してみましょう。まずはpath.plot()を使って、架空西ドイツと西ドイツのGDPの変化を確認してみましょう。

path.plot(synth.res    = Main_Fit,
          dataprep.res = Main_Data,
          Xlab = "Year", Ylab = "GDP per Capita")

これを{ggplot2}で作図する場合、以下のようなコードになります。

Synth_Result <- data.frame(
  Year = 1960:2003,
  WG   = Main_Data$Y1plot[, 1],
  C_WG = (Main_Data$Y0 %*% Main_Fit$solution.w)[, 1]
  )

Path_Plot <- Synth_Result %>%
  ggplot(aes(x = Year)) +
  geom_vline(xintercept = 1990, linetype = 2) +
  geom_line(aes(y = WG, linetype = "西ドイツ"), size = 1) +
  geom_line(aes(y = C_WG, linetype = "架空の西ドイツ"), size = 1) +
  annotate("text", x = 1990, y = 24000, label = "ドイツ統一 → ",
           hjust = 1, family = "HiraKakuProN-W3", size = 5) +
  scale_linetype_manual(values = c("西ドイツ"       = 1,
                                   "架空の西ドイツ" = 2)) +
  labs(x = "年", y = "一人当たりGDP", 
       color = "", linetype = "") +
  theme_bw(base_family = "HiraKakuProN-W3") +
  facet_wrap(~"Abadie et al. (2015)の図2") +
  theme(legend.position = "bottom",
        text = element_text(size = 16))

Path_Plot

gaps.plot()は架空の西ドイツと実際のドイツ間のGDPの差分を示したものです。

gaps.plot(synth.res    = Main_Fit,
          dataprep.res = Main_Data,
          Xlab = "Year", Ylab = "Gap in GDP per Capita")

{ggplot2}を使うと以下のようなコードになります。

Gaps_Plot <- Synth_Result %>%
  mutate(Gap = WG - C_WG) %>%
  ggplot(aes(x = Year, y = Gap)) +
  geom_vline(xintercept = 1990, linetype = 2) +
  geom_hline(yintercept = 0,    linetype = 2) +
  geom_line(size = 1) +
  annotate("text", x = 1990, y = -1000, label = "ドイツ統一 → ",
           hjust = 1, family = "HiraKakuProN-W3", size = 5) +
  labs(x = "年", y = "GDPギャップ", 
       color = "", linetype = "") +
  coord_cartesian(ylim = c(-4000, 1000)) +
  theme_bw(base_family = "HiraKakuProN-W3") +
  facet_wrap(~"Abadie et al. (2015)の図3") +
  theme(legend.position = "bottom",
        text = element_text(size = 16))

Gaps_Plot

CausalImpact

最後にCausalImpactの実装方法について解説します。必要なパッケージは{CaucalImpact}のみです。事前にinstall.packages("CausalImpact")を実行し、パッケージをインストールしてください。

{CausalImpact}パッケージの使い方

まずは{CausalImpact}パッケージと実習用データを読み込みます。実習用データは授業のサポートページ、または授業用Slackからダウンロードできます。

library(CausalImpact)

# 実習用データの読み込み
Stock_Data <- read_csv("Data/Day3_Data4.csv")

データの中身を確認します。{DT}パッケージのdatatable()関数を使用します。

datatable(Stock_Data)
変数名 説明
Date 日付
Time 時間 (一分単位)
Trend トレンド変数。13時を1とし、1分毎に1増えます
Stock_J 株価（日本）
Stock_K 株価（韓国）

NHKから正確な報道時間は不明ですが、14時6分頃だったような気がします。データから14時6分のTrend値を確認しておいてください。まずは、株価の推移を可視化してみましょう。横軸はTrend、縦軸はStock_Jにした折れ線グラフを作成します。また、14時6分（Treat == 67）に垂直線を追加すると、NHK報道の前後比較がしやすくなります。

Stock_Data %>%
  ggplot() +
  geom_vline(aes(xintercept = 67, 
                 linetype = 'NHK report of the resignation')) +
  geom_line(aes(x = Trend, y = Stock_J), size = 1) +
  labs(x = "Trend (0 = 13:00:00)", y = "Nikkei Index (Aug. 28, 2020)", 
       color = "", linetype = "") +
  scale_linetype_manual(values = c('NHK report of the resignation' = 2)) +
  theme_light() +
  theme(legend.position = "bottom",
        text            = element_text(size = 16))

報道直後において株価が急落したことが分かります。それでは、CausalImpact推定を行います。使う関数はCausalImpact()関数であり、必須引数はdata、pre.period、post.periodの3つです。他にも引数はいくつかありますが、詳細は?CausalImpactから確認してください。

data引数は第一引数であり、ベクトル、またはデータフレーム（tibble）型を指定します。こちらは株価のデータを指定します。data = Stock_Data$Stock_Jでも良いですし、Stock_DataからStock_J列のみ抽出し、これをパイプ（%>%）でCausalImpact()へ渡すことも可能です。ここでは後者の方法を採用します。

続いて、pre.periodとpost.period引数ですが、こちらは長さ2の数値型ベクトルを指定します。pre.period = c(1, 67)ならdataの1番目から67番目の要素は処置前の期間であることを意味し、post.period = c(68, 121)は68番目から121番目の要素が処置後の期間であることを意味します。

CI_wo_Cov <- Stock_Data %>%
  dplyr::select(Stock_J) %>%
  CausalImpact(pre.period  = c(1, 67), 
               post.period = c(68, 121))

結果を確認するにはsummary()関数を使います。

summary(CI_wo_Cov)
## Posterior inference {CausalImpact}
## 
##                          Average          Cumulative        
## Actual                   22881            1235580           
## Prediction (s.d.)        23307 (2.7)      1258599 (146.2)   
## 95% CI                   [23302, 23313]   [1258311, 1258875]
##                                                             
## Absolute effect (s.d.)   -426 (2.7)       -23019 (146.2)    
## 95% CI                   [-431, -421]     [-23295, -22731]  
##                                                             
## Relative effect (s.d.)   -1.8% (0.012%)   -1.8% (0.012%)    
## 95% CI                   [-1.9%, -1.8%]   [-1.9%, -1.8%]    
## 
## Posterior tail-area probability p:   0.00102
## Posterior prob. of a causal effect:  99.89817%
## 
## For more details, type: summary(impact, "report")

Absolute effect行のAverage列を見ると-426 (2.7)という結果が得られています。これはNHKの報道によって株価が約426ポイント下がったことを意味します。これは処置後期間における差分の平均値となります。Cumulateve列は差分の累積となります。他にもRelative effect列から株価が何%下がったらかが確認できます。

これらの結果を可視化するにはplot()関数を使います。

plot(CI_wo_Cov)

Originalプロットの実践は観測値、破線は予測値とその信用区間となります。13時から14時ころまでの株価は非常に安定しており、もしNHKがなかった場合の株価、つまり14時8分以降の予測値もほぼ推定になっています。そして、Pointwiseプロットは観測値と予測値の差分を示したグラフです。最後にCumulativeはPointwiseの累積です。

共変量の追加

つづいて、共変量を追加する簡単な例を紹介します。ここでは同じ時間帯を共有し、株価の変動がかなり似ている韓国の株価を使います。まずはトレンドを確認してみましょう。韓国の株価はStock_K列に入っています。

Stock_Data %>%
  ggplot() +
  geom_vline(aes(xintercept = 67, 
                 linetype = 'NHK report of the resignation')) +
  geom_line(aes(x = Trend, y = Stock_J, color = "Nikkei"), size = 1) +
  geom_line(aes(x = Trend, y = Stock_K, color = "KODEX"), size = 1) +
  labs(x = "Trend (0 = 13:00:00)", y = "Stock Price", 
       color = "", linetype = "") +
  scale_linetype_manual(values = c('NHK report of the resignation' = 2)) +
  theme_light() +
  theme(legend.position = "bottom",
        text            = element_text(size = 16))

韓国も日本と同じタイミングで株価が多少下落しましたが、全体的には安定しているように見えます。それでは、韓国の株価を追加し、架空の株価を作成してみましょう。共変量の追加はdata引数のデータフレームにStock_K列を加えるだけです。第1列は予測の対象となり、第2列以降は共変量です。したがって、1列目はStock_J、2列目はStock_Kになるようにします。

推定が終わったら結果を確認してみましょう。

CI_w_Cov <- Stock_Data %>%
  dplyr::select(Stock_J, Stock_K) %>%
  CausalImpact(pre.period  = c(1, 67), 
               post.period = c(68, 121))

summary(CI_w_Cov)
## Posterior inference {CausalImpact}
## 
##                          Average           Cumulative        
## Actual                   22881             1235580           
## Prediction (s.d.)        23299 (1.6)       1258143 (84.9)    
## 95% CI                   [23296, 23302]    [1257970, 1258304]
##                                                              
## Absolute effect (s.d.)   -418 (1.6)        -22563 (84.9)     
## 95% CI                   [-421, -415]      [-22724, -22390]  
##                                                              
## Relative effect (s.d.)   -1.8% (0.0067%)   -1.8% (0.0067%)   
## 95% CI                   [-1.8%, -1.8%]    [-1.8%, -1.8%]    
## 
## Posterior tail-area probability p:   0.00117
## Posterior prob. of a causal effect:  99.88345%
## 
## For more details, type: summary(impact, "report")
plot(CI_w_Cov)

  1. {pacman}パッケージのp_load()はパッケージがインストールされていない場合はインストール後、自動的に読み込んでくれます。インストール済みの場合はパッケージを読み込むだけですので、非常に便利なパッケージです。↩︎

  2. Alexander Fouirnaies and Hande Mutlu-Eren. 2015. “English Bacon: Copartisan Bias in Intergovernmental Grant Allocation in England.” Journal of Politics, 77 (3): 805-817↩︎

  3. 今回の例だと、地方議会ダミー変数です。↩︎

  4. 少し前まではクラスター標準誤差を計算する際、{multiwayvcov}のcluster.vcov()関数と{lmtest}のcoeftest()関数を組み合わせて使用しました。非常に面倒でしたが、今は{estimatr}のlm_robust()一つだけで解決できるようになりました。↩︎

  5. class(tidy(Bacon_Fit1_P1)[2, ])とdim(tidy(Bacon_Fit1_P1)[2, ])で確認できます。↩︎

  6. 他にも論文では人口 (対数)、選挙年ダミーも含まれていますが、ここでは省略します。↩︎

  7. 第7〜23回参院選のデータは参議院議員通常選挙データベースからダウンロードしたものを架空したものであり、第24・25回データは宋が手入力したものです。↩︎

  8. {gsynth}でもほぼ同じ結果を再現することはできます。↩︎

  9. Abadie, Alberto, Alexis Diamond, and Jens Hainmueller. 2015. “Comparative Politics and the Synthetic Control Method.” American Journal of Political Science, 59 (2): 495–510.↩︎